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Abstract 

In this study, we first present an improved version of the classical sixth-order com¬ 
bined compact difference (CCD6) scheme to enhance the convective stability of advection 
equations through an increased dispersion accuracy. This improved fifth-order dispersion- 
relation-preserving combined compact difference scheme (DRPCCD5) has been rigorously 
analyzed through the dispersion, phase speed anisotropy and stability analyses. We then 
couple the DRPCCD5 scheme with the previous fifth-order compact-reconstruction weighted 
essentially non-oscillatory (CRWEN05) scheme using a novel hybrid strategy based on 
the monotonicity-maintenance criteria. To verify the resulting ’’optimized” hybrid scheme 
(ODRPCCD5), several benchmark problems with available exact solution are investigated. 

The comparison to the previous fifth-order WENO (WEN05) scheme shows that the ODR- 
PCCD5 avoids numerical oscillation around discontinuities, handles large gradients well, 
and is much faster at the same accuracy because a coarser mesh can be used. 
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1 Introduction 


Numerical simulations of advection equations are commonly found in many applications 
of practical importance, such as shock waves, shallow water flow, magnetohydrodynamics, 
and two-phase flow models. When numerically solving such convection-dominated partial 
differential equations (PDEs), it is desirable to minimize the indispensable dispersion error, 
which is defined as the discrepancy between the numerical and actual wavenumbers, be¬ 
cause this enhances convective stability and allows for accurate capturing of small length 
scales in the wave phase [T]. 

For this purpose, dispersion-relation-preserving (DRP) approaches have been devel¬ 
oped to enhance convective stability by rigorously preserving the dispersion relation 
Furthermore, compact difference schemes offer spectral accuracy with fewer grid points 
to improve convective stability ITT UTOI . These compact difference schemes have been ex¬ 
tended to combined compact difference schemes (CCD) ITTI . in which first and second 
derivative terms are simultaneously evaluated in an implicit manner, making the scheme 
more compact and accurate. CCD schemes suffer from stability issues of boundary condi¬ 
tions when solving the PDE. In fact, these schemes need special treatment at the boundary 
nodes, in particular when simulating thin boundary layer problems. Hence, the boundary 
closures have been improved lIT^ to obtain better numerical properties, and the correspond¬ 
ing dissipation and de-aliasing properties have been discussed ifT^ . 

High spectral resolution schemes, such as the compact difference and CCD schemes, in¬ 
evitably produce numerical oscillations near discontinuities and lead to failure of the flow 
simulation. In order to avoid numerical oscillations, high resolution schemes often use 
flux/slope limiters to bound the solution gradient around shocks or discontinuities II141I15II . 
Some representative schemes belonging to this class of methods include the essentially 
non-oscillatory (ENO) scheme II161I17II and their weighted variants, known as the weighted 
ENO (WENO) |[T^IT9l . It’s well known that the ENO and WENO schemes may be too 
dissipative for compressible turbulence simulations and aero-acoustics problems. Hence, 
the compact-reconstruction weighted essentially non-oscillatory (CRWENO) scheme If20l 
has been presented, in which compact sub-stencils are identified at each interface and com¬ 
bined using the WENO weights. WENO schemes have been intensively used for problems 
containing both shocks and complicated smooth solution structures II211I22II . 

Algorithms with high accuracy are required to capture small wavelengths and non- 
oscillatory behaviors across discontinuities like shock waves. For this purpose, special 
finite difference schemes have been introduced 1231. Also, the hybrid finite difference 
scheme based on the minimized dispersion and controllable dissipation (MDCD) technique 
has been developed to solve advection equations. This MDCD technique has been coupled 
with an optimized WENO scheme to make discontinuity capturing possible |[24l . Many 
researchers have also proposed various alternative ways to improve the numerical schemes 
Il25l - l2^ . However, accuracy still remains a challenge because, to our knowledge, most 
if not all existing numerical schemes suffer from the drawback that they switch to a non¬ 
compact scheme at and near discontinuities, resulting in a loss of resolution. 

In this study, a fifth-order dispersion-relation-preserving combined compact difference 
(DRPCCD5) scheme which has better DRP properties than previously reported compact 
difference schemes over a considerable range of wavenumbers is proposed. This scheme 
ensures that resolved energy components propagate closer to the correct physical speed, 
and that complex phenomena, involving interactions among different wavelength scales, 
can be captured. Furthermore, the DRPCCD5 scheme is coupled with the CRWEN05 
scheme using a novel hybrid strategy based on the monotonicity-maintenance criteria. The 
numerical properties of the resulting ’’optimized” hybrid scheme (ODRPCCD5) are then 
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rigorously analyzed using several benchmark problems. 

This paper is organized as follows. Section 2 describes discretization of a standard 
advection equation and the time marching method, which is used in the present study. 
The schemes construction is carried out in Section 3. Section 4 includes the fundamental 
analysis of dispersion, dissipation, phase speed anisotropy, numerical group velocity, and 
numerical phase velocity for the proposed DRPCCD5 scheme. Several benchmark tests are 
performed in Section 5 to validate the ODRPCCD5 scheme. Section 6 draws concluding 
remarks based on the results presented in Section 5. 


2 Time marching method 


The one-dimensional linear wave equation can be expressed as 

du df 

-h — = 0. 

dt dx 


( 1 ) 


where t is time, x the spatial coordinate, u the field variable, and f = cu with c the constant 
propagation speed of the wave. A conservative finite difference discretization of Eq. O 
results in an ordinary differential equation, which can be expressed as 


dui 

dt 


= Fi{u) 


h^fi+k 



( 2 ) 


where h is the grid spacing and is the numerical approximation of flux between points 
Xi and x,+i. In the present study, we apply the fourth-order Runge-Kutta (RK4) scheme 
and the sixth-order symplectic Runge-Kutta (SRK6) scheme 1291 for time evolution. The 
explicit RK4 scheme reads 


„(2) = ^,(1) + ^ +f («('))), 

„(4) =„(3) + ^ («(!))-4F(M(2))+i7(„(3)))_ (3) 

6 

For the SRK6 scheme, given the solution at t = nAt, the solution is obtained from 
the following iteration. We start with computing and = F{u^j^), where j=\ to 3, 
by numerically solving the following equations iteratively: 

»">=«" +4< [Tf(') + [l + |)fP) + 2 + P)], (4) 

=„"+A, [(2 _ 2)f(i)+ ^(2 )+(^+(3)], (5) 

„(5' = „"+A, 1(2 _ £)f(1) + ( a |)f( 2)+2f (3)], (6) 


where c = ^ y These updated values correspond to the times t = n + {^ +c)At, 1 = 

n + ^At, and t =n + — c)At, respectively. Upon reaching the user’s specified tolerance 

(10^^), the solution dlt = {n+ l)At is obtained as 


/+! + ^ [^/ 7 ( 1 )+ 4 / 7 ( 2 ) + 5 / 7 ( 3 )]. 

9 ^2 2 ^ 


(V) 
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The RK4 scheme is mainly used to run the numerical tests in this study because the implicit 
SRK6 scheme provides nearly the same results, but is very time-consuming (see results of 
linear advection problem #1 in Table 1 and Fig. 7). 


3 Numerical Schemes for spatial discretization 


3.1 Fifth-order non-compact difference scheme 

The numerical flux can be reconstructed using a left or right biased interpolation If20ll . 
The appropriate interpolation is chosen based on the sign of the wave speed, which in the 
case of a scalar PDF is given by 


^!+ 1/2 — 0 , 

= 127^1/2 if c,-+1/2 <0. (8) 

where the superscripts + and — denote left and right biased interpolations respectively. 
Note that the approximation of the left biased numerical m,+i /2 is described in this section. 
The first derivative term ^ can be approximated to the desired order (r), reading 


du 

dx 


lx=Xj 


“i+1/2 — Mi-1/2 

h 




where 12;+ 1/2 for odd r is computed using the linear reconstruction on a stencil 


(9) 


(r-l)/2 

Mi+1/2 = hUi+k- (10) 

k=-{r-\)/2 


Here bk is the coefficient and Ui+k = 

_ 1 

Mi + 1/2 — ^“‘-2 — 


Uxi+kh- For r = 5, it reads 
13 47 27 1 


60 


-Uj- 




20 


Mi+2- 


( 11 ) 


This scheme has fifth-order spatial accuracy according to the derived modified equation 
given below 


du 

dx 


du 

\ exact 
OX 




( 12 ) 


The drawback of this scheme is that the magnitude of the leading error of the resulting 
scheme is too large. In addition, this scheme must adopt a very fine mesh to correctly 
capture the important advection flow structures. A space-time accurate numerical sim¬ 
ulation of advection problems requires higher spatial resolution and dispersion-relation- 
preserving (DRP) properties. Such schemes act as an important numerical tool to solve 
complex physical problems displaying a large bandwidth of spatio-temporal scales. For 
this reason, we develop a new fifth-order dispersion-relation-preserving combined compact 
difference (DRPCCD5) scheme in the Section 3.2. 


3.2 Fifth-order dispersion-relation-preserving combined com¬ 
pact difference (DRPCCD5) scheme 

In this section, we present an improved upwind combined compact difference scheme. 
The first and the second derivative terms ( and in a four-point grid stencil are ap- 
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proximated as 


du , 


du , 


du , 


i-i + ^ + ;+i 

ax ax ax 


= J_{c\Ui-2 + C 2 Ui-i +C 2 Ui)-h[ +^ 233 l( + ^333 




d\ 


d^u , 


dx^ 
du, 


’ dx^ 


du , 


|(+i 


Id^u d^u Id^u. 3 . , 9 I 


( 13 ) 


(14) 


The coefficients shown in Eq. (flTh are derived through Taylor series expansion. Elimi¬ 
nation of the leading truncation error terms in the modified equation analysis enables us to 
get the formal accuracy order of six ifTll . 

Derivation of the coefficients in Eq. (fT3]) is started from performing Taylor series ex¬ 
pansion on the terms m,_ 2 , m/-i, ^\i~u ^\i, ^\i+u ^\i and with respect 

to Ui to get the modified equation. The six leading truncation error terms derived in the 
modified equation analysis are then eliminated to get a set of six algebraic equations 


Cl -hC2-hC3 
—2ci — C2 — a\ — 

Ac\ -|-C 2 + 2ai — la^, —2b\— 2b2 — 2bi, 
8ci -|- C2 + 3ai -|- 3(33 — 6t7i -|- Gb^, 

16ci -|-C 2 -|-4ai — 4a3 — I2bi — \2b2 
32ci 4“ C2 “h 5cti 4“ 5 ce3 — 20b[ -4 20Z?3 


= 0, 

(15) 

= 1, 

(16) 

= 0, 

(17) 

= 0, 

(18) 

= 0, 

(19) 

= 0. 

(20) 


Derivation of two further algebraic equations are needed to determine all eight coefficients 
in Eq. (fTSl) . One way of deriving the two equations so as to get a better approximation of 
is to reduce numerical error of the accumulative type. We can then expect to retain the 
theoretical dispersive property of ^ l!^ . 

Our strategy of achieving the goal of reducing numerical dispersion error is to match the 
exact and numerical wavenumbers. Use of this underlying approach amounts to equating 
the effective wavenumbers a and a to those shown on the right-hand sides of Eqs. (|2T1) 
and (1221) {21. Eollowing this line of derivation, we are led to get the two equations for a'h 
and a"h as follows 

iah + 1 +^ 3 ^“) = (cjc^^m + ^ 26 -*“^' + C 3 ) - (ia"/i)2()qc-‘“'’ +^72 + 

( 21 ) 


{iahf{--e-^^ + 1 - -c“) = (3c-‘“'' - 6 + 3c“) - iah (22) 

8 8 8 8 

Equations (|2T1) and (1221) are solved to get the expression for a h which has been used sub¬ 
sequently to minimize the dispersion error. The real and imaginary parts of a h provide 
information regarding the dispersion error (phase error) and dissipation error (amplitude 
error), respectively. 

To improve the dispersive accuracy for a', the exact value aJi should be very close to 
91[a'/7], where 91[a'/i] denotes the real part of a'h. To achieve the goal of improving solution 


5 


accuracy, the positive-value error function E{a) defined below should be very small over 
the following integration interval for the modified wavenumber ah 


E(a)= J ' [W-{ah-^i[a'h])Yd{ah). (23) 

In Eq. (|2^ fhe weighting function W is chosen fo be fhe denominator of {ah — Si[a'h]). 
This choice facilitates us to integrate E (a) exactly. To make the error function defined in 
0 < a/j < ^ fo be positive and minimal, fwo exfreme conditions given by 



(24) 



(25) 


are enforced. These fwo consfrainf equations enforced for maximizing fhe dispersion ac¬ 
curacy are used fogefher wifh fhe ofher six algebraic equafions derived from fhe modified 
equation analysis fo gef nof only a smaller dissipation error buf also an improved dispersion 
accuracy. Nofe fhaf several infegrafion ranges have been numerically defermined so as fo 
find fhe besf one fhaf renders fhe smallesf value of E. 

The resulfing eighf infroduced unknown coefficienls can be defermined as a\ = 0.8873686, 
fl 3 = 0.0491178, bi = 0.1495320, b2 = -0.2507682, b^ = -0.0123598, ci = 0.0163964, 

C 2 = —1.9692791 and C 3 = 1.9528828 from fhe above reduction of dispersion and dissipa¬ 
tion errors. The upwinding scheme developed fheorefically in four stencil poinfs / — 2, / — 1, 
i and i -I- 1 for has fhe spatial accuracy of order fiffh according fo fhe derived modified 
equafion given below 


^ + 0.0000077381655315119445 h^^+H.O.T. . (26) 

(JJv (JJv yJJv 

If is nofed fhaf, unlike our sfrafegy, Zhou el al. |[28]l chose fhe coefficienl C 3 as free parame- 
ler so fhaf fhe ofher seven coefficienls are expressed as fhe linear funcfions of C 3 by Taylor’s 
expansion. Then Ihese eighf coefficienls were numerically defermined by fhe sfandard se¬ 
quential quadratic programming (SQP) melhod |[30l fo minimize fhe error function shown 
in Eq. (1231) . However, Ihis oplimizafion resulf is highly sensitive fo fhe inifial guess of C 3 , 
as pointed by Zhou el al. If28]l . 

Define firsl fhe values of u af fhe half nodal poinfs / ± ^ as follows: 


“r+l/2 — l\Ui-\ +Y2te— [(0 CiM,-1/2 + CX2M,'+3/2) +h{^l^i-l/2 + P2“i+l/2 + 

(27) 


and 

Ui_i/2 = YiWi-2+Y2te-1 — [(ctlte-3/2+ Ct2te+l/2) +^(Pl“,--3/2 + p2“i-l/2 + p3“!+l/2)]- 

(28) 


One can Ihen subslilule Ihem info Eq. (|9ll fo gef fhe algebraic equafion for al fhe node i. 
Derivalion of OC/, P, and y, is fhen followed by comparing fhe coefficienls derived in Eq. (fT3]) 
for ^|;. Affer a lerm-by-lerm comparison of Eq. I®, we are led fo gel fhe coefficienls 
as follows: oci = 0.8873686, a 2 = 0.0491178, Pi = 0.1495320, P 2 = -0.2507682, P 3 = 
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-0.0123598, Yi = -0.0163964, \ = 1.9528828. In brief, of DRPCCD5 scheme for 
Q+ 1/2 >0is given by 

= -0.0163964m,•_! + 1.9528828/,- - [(0.8873686m,-_i/2 +0.0491 178m,-+3/2) 

+/j(0.1495320m4i/2 -0.2507682m4i/2 -0.0123598m43/2)]- 

(29) 

Thus, the magnitude of the leading error term in the compact interpolation is less than the 
corresponding non-compact interpolation on the same order (see Eqs. (fT^ and (1261) ). m,+i /2 
of DRPCCD5 scheme for c,-+i /2 < 0 can be similarly derived: 

= 1.9528828m,-_i -0.0163964m,- - [(0.0491 178m,-_ i/2 +0.8873686 m,-+3/2) 

+/j(0.0123598m4i/2 +0.2507682m,.+i/2 -0-1495320m,.+3/2)]. 

(30) 


3.3 Weighted essentially non-oscillatory (WENO) scheme 

Advection equations admit discontinuous solutions. Weighted essentially non-oscillatory 
schemes are designed to achieve the high order of accuracy at smooth regions and switch 
to lower order interpolation to avoid oscillations near discontinuities. 

3.3.1 Fifth-order WENO (WEN05) scheme 

The form of the interface flux reconstructed by the WEN05 scheme 1191 reads 

fj+l /2 = g (7(Ol + 02)//'-1 + -(llCOl + 5(02 + 2 ( 03 )/;- 

+ -(2(02 + 5(03)/;-+!-^/;+2- (31) 

6 6 

In the above equation, we write 

OCf ~ Ck 

Oik = ^, o.u= ~ , ^ = 1,2,3. (32) 

l.k^k (pyt + e)^ 

The optimal weights are ci = ^, C 2 = ^ and 03 = ^. A very small number (e = 10^^) 
is used to prevent division by zero. The smoothness indicators ^k are given to detect large 
discontinuities and automatically switch to the stencil that generates the least oscillatory 
reconstruction by 

Pi = {^(/-2 -2/--1 +/-)V ^(/--2 -4/--1 +3/-)2, 

P 2 = {^(/-i -2/- + /-+i)H^(/--i -/-+i)", (33) 

P 3 = Y 2 “ 2/-+1 + fi+if + ^ (3/- — 4/,-+i + 

The WEN05 scheme gives fifth-order accurate results in smooth regions of the solution 
and is non-oscillatory near discontinuities. 
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3.3.2 Fifth-order compact-reconstruction WENO (CRWEN05) scheme 

The drawback of higher order WENO schemes is the increasingly wide stencil when 
increasing the order of accuracy. Therefore, the CRWEN05 has been constructed using 
three third-order compact interpolations as candidates GOl . The CRWEN05 scheme can 
be expressed as 


21 /^ 12 l/^ 

kj®! + + [ 3 ®! + 2(®2 + ®3)].^’+l/2 + 2®3.^-+3/2 

(0i 5((0i+(02) + (03 . , (02 + 5(03 . 

“ +- 6 -+- 6 - 


(34) 


Note that fi^i /2 io E(l- (1341) is the approximation of the left biased numerical flux 

for / (m)|ji:=x ,_|_,/2 ^ 0. Since the weights (Ok in Section 3.3.1 are overly dissipative, they are 
determined, as suggested in the literature II311I32II . using at as 


ak = Ckil-\ -^), k = 1,2,3. 

£ + Pi: 


(35) 


Here, the optimal weights are ci = j, C 2 = 5 and C 3 = ^. x is simply defined as the absolute 
difference between Po and P 2 . 


3.4 Fifth-order optimized dispersion-relation-preserving com¬ 
bined compact difference scheme (ODRPCCD5) 

In this section, we briefly present the hybrid strategy to couple CCD with WENO 
schemes proposed by II241I271I28I and our novel hybrid strategy based on the monotonicity- 
maintenance criteria. Both strategies are compared with each other in Section 5.1.2. 

3.4.1 Hybrid strategy by If24l[27l[28]l 

Eollow the hybrid strategy of II241I27[|2^ . the numerical flux f 1 + 1/2 can be written as 

f — fDRPCCDS , / 1 3 fCRWENOS 

/r+1/2 — (^1+1/2//+1/2 + (t “Ct(+l/2j//+i/2 ■ (3 d) 

In the above, ct,+i /2 is the weight function and its detailed formulation can be expressed as 

a,+i/2 = min(l,^), (37) 

fc 

where rc is constant and r,+i /2 is a smoothness indicator, determined as 

^(■+ 1/2 =min(n-,ri+i), 

with 

|2A/;'+i/2A/;'i/ 2|+ei 
where A/;+i /2 = /;+i -/; and 81 = lO^*’. 


(38) 


(39) 
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3.4.2 Present hybrid strategy 


We first define the monotonic range in our present hybrid strategy. The field variable 
u{x,t) is normalized by 


u{x,t) 


u{x,t) — u1^ 


(40) 


When substituting the node values u-Lj and into Eq. (l40l) . these values can be nor¬ 
malized as u'l_^ = 0 and = 1, respectively (see Fig. 1). As shown in Fig.l, we then 
establish our hybrid strategy based on monotonicity-maintenance criteria by requiring face 
values Ui+ii 2 - 


and m,_i/ 2 : 


< Mi+l/2 < 1) 

0 < 1 /2 < • 


(41) 


(42) 


where M;±i /2 is calculated by substituting face value ^,± 1/2 irito Eq. (l40l) . 

In addition, the new Ui value must be constrained to maintain monotonicity by the fol¬ 
lowing formulation 




(43) 


We discretize Eq. O as 


=u’l- Cr(M;+i /2 - M;- 1 / 2 ), (44) 

where Cr = ^. Substituting Eq. (1441) into left-hand inequality of Eq. (|4^ leads to 

“r+I/2 < “i-1/2 + ^ (“f “ ) • (45) 

Since m;_i /2 > 0 and < 0, the worst-case condition in Eq. (1451) is given by m;_i /2 = 0 
and = 0. It means that Eq. (1451) can be rewritten as 

u" 

m;+i/2 <^- (46) 

Thus, the monotonic range can be determined by Eq. (BTl) . Eq. (|4^ and 0 < ii” < 1, 
as shown in the shadow region in Fig. 2. The slope of the Courant-number-dependent 
boundary (dashed line in Fig. 2), changes with Cr. 

Once the monotonic range is defined, we then calculate and substitute it into 

Eq. (l40l) to get if,+ 1 / 2 , and estimate whether m;+i /2 locates in the monotonic range. If yes, 
set m,+i /2 = If not, set m,+i /2 = or m,+i /2 = nf. For clarity, the steps are 

given as follows: 


Step 1: if c > 0, set 

Step 7 according to Fig. 3(a). 

Step 2: If c < 0, set 

Step 7 according to Fig. 3(b). 


_ ',DRPCCD+ 
~ ^i+l/2 

_ /^.DRPCCD- 
r+1/2 


r~,CRWENO _ r,CRWENO+ 
“i+1/2 “ “;+i/2 

aCRWENO _ ,'^CRWENO- 

“i+1/2 “ “i+1/2 


and perform Step 3 to 


and perform Step 3 to 


Step 3: Compute B = «£> — mj/; if |B| < 10 set m,+i /2 = uc- 

Step 4: If |B| > 10^^, compute % = {uc — mc/)/I3; if this is less than 0 or greater than 1, 
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again set m,+i /2 = “c- 

Step 5: Compute m,+i /2 = and 

Step 6: If m;+i /2 < %, set m;+i /2 = 

Step?: If m;+i/2 > Mc/Cr, set M,-+i/2=Mc/Cr; if m,+i/2 > C reset m,+i/2 = 1- Construct 
W;+l/2 = + Uu- 

Step 8: Calculate face values ii±i /2 = C;±i/ 2 M,+i /2 and update into the next time step ac¬ 
cording to Eq. (O- 

Coupling DRPCCD5 with CRWEN05 using this hybrid strategy leads to the ODR- 
PCCD5 scheme. 


4 Fundamental analysis 

4.1 Dispersion and dissipation errors 

The solution for the model equation 


Mf + C Mx = 0, 


is given by 


u = 




(47) 


(48) 


where i = and is the Eourier mode of the wave number a. Differentiation of the 
above equation leads to 


^\exact=io.h'^e'^. (49) 

a;c h 

The approximated derivative term ^ can be similarly written as 

^ A. A. 

\numerical — lOC n ^6 — (A:^ T lA,j ^ 6 . (50) 

Here, and A„ denoting the real and imaginary parts of ah (cf. Eq. (I2TI) '). account for the 
dispersion and dissipation errors, respectively. 

Eig. 4 shows the dispersion and dissipation characteristics of fifth-order non-compact 
finite difference scheme (EDS) |[20l . fifth-order compact difference (CDS) scheme EOl . 
eighth-order optimized compact difference (OCD8) scheme |[33l and our proposed DR- 
PCCDS scheme. It can be seen that the DRPCCDS scheme has a better spectral resolution 
than the OCD8 scheme. The dispersion property of the DRPCCDS scheme is better than 
those of the other schemes because of the improved dispersive accuracy. Eurthermore, at 
frequencies with low dispersion error, the DRPCCDS scheme has less dissipation than the 
other schemes. 


4.2 Assessment of the phase speed anisotropy 

In anisotropic two-dimensional problems, first-order differencing schemes tend to pro¬ 
duce phase space errors iTTlfm . To evaluate the phase space error of our DRPCCD5 scheme, 
we take the following two-dimensional advection equation into consideration 

U,+CxUx + CyUy=0. (51) 
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Here, Cx = c cos(9) and Cy = c sin(9) denote the velocity components along the x and y 
directions, respectively. For a wave propagating at the angle 9 (= tan^^(^)) with respect 
to the A:-axis, the numerical phase speed anisotropy can be derived as follows iTTl 

cos{Q)3i[a h{ah cos{Q))]+sin{Q)3{[a h{ah sin{Q))] 
c^~ ■ 

One can find from Fig. 5 that our proposed scheme reproduces phase speed anisotropies 
much better than the sixth-order combined compact difference (CCD6) scheme ifTTTl at all 
scaled wavenumbers. 


4.3 Amplification factor, numerical group velocity and numeri¬ 
cal phase velocity 

The properties, such as amplification factor, numerical group velocity and numerical 
phase velocity fOl, of the present DRPCCD5 scheme are analyzed by solving the one¬ 
dimensional wave equation, where the fourth-order accuracy Runge-Kutta (RK4) scheme 
is applied in time evolution. The present scheme is compared with the previous sixth-order 
combined compact difference (CCD6) scheme ifTH . The general numerical solution of Eq. 
(STJl is identified as 


u(x,„,r”) = J f/(a,r")e‘“'” da, 

(53) 

such fhaf fhe inilial solulion is given by 


u{xm,t = 0) = J Ao(a)e‘“"' da. 

(54) 

Nole fhaf fhe u{xm,t") can be oblained by subsliluling fhe above inilial condilion as [91 

u{x,n,f) = 1 Ao(a)(G2 + G2)5e‘(“'"“"P) da. 

(55) 

In Eg. (1551). fhe numerical amplification factor G{a) is defined as G{a) 

— (7 1 if! — 

-Gr + lGi- 

The ferm P is oblained as lanP = — The numerical group speed 

velocity are oblained as 

and numerical phase 

Eg (a) 1 dp 

(56) 

c hCrda' 

Vp{a) P 

c coAt ’ 

(57) 


where Cr = ^ = ^ denofes fhe Couranf number. 

In Figs. 6(a) and (b), fhe amplificalion facfors are nafurally sfable over a large range of 
coA? for bofh DRPCCD5 and CCD6 schemes. Figs. 6(c) and (d) show fhe comparison of 
fhe variations of ^ in fhe ah — (oAt plane for fhe fwo numerical schemes discussed above. 
If one defines fhe area bounded by fhe confour lines of ^ = 0.95 and ^ = 1.05 as a DRP 
region, fhe DRPCCD5 scheme can resolve fhe DRP region up fo ah = 2.5, while fhe CCD6 
scheme only reaches ah = 1.68. If can be clearly seen fhaf fhe DRPCCD5 scheme has fhe 
heifer DRP properly. Figs. 6(e) and (f) give fhe conlours of fhe numerical phase speed. 
Similarly, defining fhe DRP region as bounded by ^ = 0.95 and ^ = 1.05, one can see 
fhaf fhe DRPCCD5 scheme resolves a 10% larger DRP region lhan fhe CCD6 scheme. 
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5 Numerical results 


5.1 One-dimensional problems 

The ODRPCCD5 scheme is tested to solve three linear advection equations and one invis- 
cid Burgers’ equation. The L 2 -errors and their corresponding spatial rates of convergence 
are tested for the linear advection problem#!. The computational costs are compared using 
different spatial discretization schemes for the linear advection problems#2. Two hybrid 
strategies descried in section 3.4 are used to solve linear advection problem#2. Finally, we 
extend ODRPCCD5 scheme to solve the one-dimensional Euler equations of the polytropic 
gas dynamics. 

5.1.1 Linear advection problem #1 

The problem with the smooth initial condition u{x,0) = sin(27ix) for Eq. ([Hi with c = 1 is 
solved. Periodic boundary conditions are applied at two boundaries of the region 0 < x < 1. 

To compare the computational efficiency of time evolution, we solve this problem by us¬ 
ing the sixth-order implicit symplectic Runge-Kutta scheme (SRK6) and the fourth-order 
explicit Runge-Kutta scheme (RK4). The twin-tridiagonal coefficient matrix for the DR- 
PCCD5 scheme is solved by the computationally effective solver, including twin-forward 
elimination and twin-backward substitution techniques, which is described in ifTH . All the 
computational times are obtained using a Core i7, 3.40 GHz computer with 64.0 GB of 
RAM. 

Table 1 shows that the SRK6 scheme costs more CPU time than the RK4 scheme when 
the same spatial scheme and grid are used. Pig. 7 shows that the computational errors 
mainly come from the spatial discretization, by comparing RK4AVEN05 and RK4/DRPCCD5. 
Therefore, we employ the RK4 scheme for time evolution in the following numerical cases. 
The L 2 -eiTors and their corresponding spatial rates of convergence, by using DRPCCD5, 
WEN05, and CRWEN05 schemes, are given in Table 2 with time step At = 1 x 10^^. It 
can be seen that all schemes can approximately achieve their theoretical order of accuracy. 


5.1.2 Linear advection problem #2 

The one-dimensional linear equation Ut + Ux = 0 is solved considering the following 
initial condition Il34l : 


u{x,0) 


^{G{x,z-?>) + G{x,z + 5)+4G{x,z)) 

1 

< l-|10(x-0.1)| 

g(F(x,a —5) -|-5) +4F{x,a)) 

. 0 


; — 0.8 < X < — 0.6 

; -0.4<x<-0.2 
; 0<x<0.2 

; 0.4<x<0.6 
; otherwise. 


(58) 


where G{x,z) = , F{x,a) = (max(l — a^(x —a)^,0))^'^^. The constants are taken 

as a = 0.5, z = —0.7, 5 = 0.005, a = 10, and p=(log 2)/365^. This initial condition consists 
of a discontinuous square wave, an exponential wave, a triangular wave, and a parabolic 
wave. Periodic boundary conditions are imposed here. The time step is chosen as At = 
0.05#. Pig. 8 shows the exact waveform and the waveform obtained by the WEN05 and 
ODRPCCD5 scheme on a grid with 200 points at t = 2 and t = 4. Pigs. 9 and 10 show the 
magnified solution for the exponential and square waves at t = 4. In Pig. 9, one can see 
that the ODRPCCD5 show less clipping at the extreme than the WEN05 in the case of the 
exponential wave. In Pig. 10, the ODRPCCD5 scheme is less dissipative than the WEN05 
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scheme across the discontinuities. Since the DRPCCD5 scheme is not classified to be a 
non-oscillatory scheme, the predicted kinks near the root of square wave is computation¬ 
ally inevitable. Comparing the magnitude of errors produced by WEN05, DRPCCD5 and 
ODRPCCD5 for this test problem shows that ODRPCCD5 performs better. 

The computational costs using WEN05, DRPCCD5 and ODRPCCD5 schemes are 
compared based on different grids, as shown in Table 3. The ODRPCCD5 scheme needs 
more CPU time than the other two schemes if the same grid is used because this scheme 
is hybrid. However, the spectral properties of the ODRPCCD5 scheme imply that it may 
apply a coarser grid to achieve the same resolution as the WEN05 scheme at the same order 
of convergence. As shown in Pig. 11, the ODRPCCD5 scheme with 600 grids reaches a 
better resolution than WEN05 scheme with 1600 grids. Meanwhile, it only needs 231s in 
comparison with 3.9^ by WEN05 scheme. 

The two hybrid strategies introduced in section 3.4 are used to solve the advection 
equations, and the numerical results are plotted in Pig. 12 with 400 grids and At = 0.05h 
at t = 2.0. In Pig. 12, we can see that the solution is not damped when using the previous 
hybrid strategy by Il24ll27ll2^ when = 0.1. Therefore, this hybrid strategy needs an 
appropriate trial parameter (r^) to damp the oscillation. In contrast, our hybrid strategy 
is based on the monotonicity-maintenance criteria, which automatically limits oscillations 
and captures discontinuities, as shown in Pig. 12(b). 

5.1.3 Linear advection problem #3 

We solve the linear equation Ut + Ux = 0, —1 < x < 1, with periodic boundary condition 
If35]l . The initial condition reads 

( ; —l<x<—^ 

u{x,t=0) = <. |sin(27tx)| ; -i<x<0 (59) 

[ 2x—1 — gsin( 37 tx) ; ^<x<l 

The predicted results in the domain with 200 grid points are plotted in Pig. 13 at t = 20. It 
can be seen that ODRPCCD5 scheme performs better than the WEN05 scheme. 

5.1.4 Non-linear advection problem 

We solve the Burgers’ equation n, -|- {0.5u^)x = 0, — 1 < x < 1, with periodic boundary 
condition. The initial condition is m(x,0) = 2-|-sin(7t(x-|- 1)). The solution to Burgers’ 
equation is smooth for t < ^ and it develops shocks for t = ^ . The results obtained at 
t = 0.3 (before shock) and t = 0.35 (after shock) are plotted in Pig. 14 in the domain with 
200 grid points. The time step is chosen as At = O.lh in this computation. We observe that 
ODRPCCD5 gives better results than the DRPCCD5 scheme at f = 0.35. 


5.1.5 The Shu-Osher problem 

In this case, we solve the one-dimensional Euler equations of gas dynamics |[24]| 



p = {y-l){E--pq 


pq \ 


pq^+p =0, 

(60) 

.^(E + p)/ 


\pq^), Y= 1-4. 

(61) 
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where p, q, p and E are the density, velocity, pressure and total energy of the conserved 
fluid, respectively. The initial conditions are 


{p,u,p) = 


(3.857143,2.629369,10.3333) 
(1 +0.2sin(5x),0,1) 


if X <1 
otherwise 


(62) 


This test case leads to very strong shock waves and is employed to validate the shock¬ 
capturing capability of the proposed ODRPCCD5 scheme. Reflective boundary conditions 
are applied at both x = 0 and a: = 10. Since the exact solution for this problem is not 
available, the solution computed in 10000 grids is considered as the exact solution. Fig. 15 
shows waveforms at t = 0.45, t = 0.9, t = 1.35 and f = 1.8 (grid spacing h = ■^, time 
step At = 0.05/j). It can be seen that the shock-waves are well reproduced by our proposed 
ODRPCCD5 scheme. 

5.2 Two-dimensional problems 

In this subsection, we illustrate the capacity of the ODRPCCD5 scheme through two- 
dimensional numerical simulations. 

5.2.1 Vortex flow problem 

The equation (|)f -|- {u(^)x + (v(|))v = 0 is solved using an initial circle shape in a square of 
unit length, within which the vortex flow field («, v) is given by If36]l 


= —sin^(7tx)sin(27ty), 
V = sin^(7ty)sin(27tx). 


(63) 

(64) 


u = 


The radius of the circle is 0.15 located at the center (0.5,0.75). At t = T the flow field was 
reversed, so that the exact solution att = 2T should coincide with the initial condition. This 
problem has been known to be computationally challenging since its solution is stretched 
and tom by the vortex flow where a very thin filament having a scale of single mesh size 
can be generated. 

Computations were performed for T = 2.5 and At = The predicted results of 
WEN05 and ODRPCCD5 are compared for the calculation of (|) = 0. The results obtained 
in 100 X 100 grids at t = 1.5, 2.5, 4,5 are plotted in Fig. 16. It is clear that the solution 
computed using the ODRPCCD5 scheme is maintained within a thin and elongated fila¬ 
ment on the scale of one grid spacing. On the contrary, the WEN05 scheme results in a 
considerable reduction of the area at the head and tail of the filament. In Fig. 16(d), one 
can see that the solution computed using our proposed scheme returns to its initial state. In 
Fig. 17, the ODRPCCD5 scheme using a 100 x 100 mesh can reach the same resolution 
att = 2.5 as the WEN05 scheme using a 200 x 200 mesh. Hence, the ODRPCCD5 needs 
less CPU time (8.30^) than the WEN05 scheme (17.325). 

5.2.2 Zalesak’s problem 

The Zalesak’s problem II371I38II is one of the best known benchmark cases for testing the 
developed advection scheme. The slotted disk has a radius of 15 and a slot width of 5. It is 
initially located at (50,75) in the domain of size (100,100). The prescribed velocity field is 
given as 


^ _7l(50-y) _7t(x-50) 

“ 314 “ 314 


(65) 
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The results predicted for 100 x 100 grid points at t = 5071, t = IOOti, t = 15071 and t = 
20071 are plotted in Fig. 18(a). The results are also plotted in Fig. 18(b) in the domain 
with 200 X 200 grid points. The solution computed with the proposed scheme is in good 
agreement with the exact (or initial) solution as shown in Fig. 18(b). 


6 Concluding remarks 

In this paper, a fifth-order dispersion-relation-preserving combined compact difference 
(DRPCCD5) scheme has been proposed, which shows increased dispersion accuracy and 
improved dispersion-relation-preserving properties compared to the CCD6 iTll scheme. To 
make discontinuity capturing possible and handle large gradients, an optimized DRPCCD5 
scheme (ODRPCCD5), which couples the DRPCCD5 and CRWEN05 schemes, is con¬ 
structed using a novel hybrid strategy based on the monotonicity-maintenance criteria. The 
numerical solutions of linear problems show that our ODPRCCD5 scheme performs very 
well and is much faster than the previous WEN05 scheme at the same accuracy. In ad¬ 
dition, the ODPRCCD5 scheme produces non-oscillatory solutions of the Euler equations 
in domains with discontinuities, and it can handle sharp resolutions when solving the two- 
dimensional vortex flow and Zalesak’s problems. We plan to apply our algorithm to solve 
the three-dimensional Navier-Stokes equations for the simulation of two-phase flows in 
future studies. 
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Scheme 

grids 

CPU times ( 5 ) 

SRK6/DRPCCD5 

40 

22.93 


80 

42.13 


160 

79.03 


320 

145.29 

RK4/DRPCCD5 

40 

8.90 


80 

13.11 


160 

22.40 


320 

41,46 

RK4AVEN05 

40 

5.83 


80 

6.95 


160 

9.68 


320 

15.39 


Table 1: Comparisons of the computational costs for the different schemes at t = 1000 with 
At = 0.001. This problem is described in section 5.1.1. 


Scheme 

grids 

L 2 error norms 

rates of convergence 

WEN05 

20 

3.724x10-^ 



40 

8.297x10-6 

5.488 


60 

8.832x10-^ 

5.524 


80 

1.835x10-^ 

5.461 

CRWEN05 

20 

8.056x10-6 



40 

1.198x10-^ 

6.071 


60 

1.229x10-^ 

5.615 


80 

2.501 xl0-‘^ 

5.534 

DRPCCD5 

20 

1.207x10-6 



40 

2.597x10-^ 

5.539 


60 

2.783xl0-‘^ 

5.508 


80 

5.750x10-^0 

5.482 


Table 2: The predicted L 2 -error norms and the corresponding spatial rates of convergence for 
the solutions predicted with SRK6 scheme at t = 1 in a domain containing four chosen meshes. 
This problem is described in section 5.1.1. 
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Scheme 

grids 

CPU times ( 5 ) 

WEN05 

200 

0.078 


400 

0.32 


600 

0.59 


800 

1.04 


1600 

3.91 

DRPCCD5 

200 

0.20 


400 

0.79 


600 

1.68 


800 

2.99 


1600 

11.43 

ODRPCCD5 

200 

0.26 


400 

1.21 


600 

2.37 


800 

4.99 


1600 

19.20 


Table 3: Comparisons of the computational costs for the different schemes at t = 4 with Courant 
number 0.05. This problem is described in section 5.1.2. 


20 




(a) 

Figure 2: Monotonic range and normalized variable values. The dashed line is a Courant- 
number-dependent slope of 
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(b) 


Figure 3: Definition of upstream (U), downstream (D) and central (C) node-values, (a) c > 0; 
(b) c < 0. 
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(a) 



(b) 

Figure 4: Comparison of Ki{ah) and Kr{(xh) amongst the proposed DRPCCD5 scheme, 
WEN05 scheme lIT^ . CDS scheme [l20ll . and OCD8 scheme lf3^ . (a) Kt; (b) Kr. 
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(b) 

Figure 5: Comparison of the predicted phase speed anisotropy, which is plotted against 0, for 
the proposed DRPCCD5 scheme and the CCD6 scheme of Chu and Fan ffTTl . (a) DRPCCD5 
scheme; (b) CCD6 scheme ifTTll . 
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(a) 


(b) 




(c) (d) 




(e) (f) 

Figure 6: Amplification factor (a) and (b), scaled numerical group speed (c) and (d), and scaled 
numerical phase velocity (e) and (f) contours for RK4 time-integration scheme with: (a)(c)(e) 
present DRPCCD5 and (b)(d)(f) CCD6 scheme IfTTI . 
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(a) 

Figure 7: The predicted results for linear advection problem#! are plotted using 40 grids at 
t = 1000. 
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WEN05 
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(a) 


WEN05 

0DRPCCD5 



(b) 

Figure 8: The predieted results for linear adveetion problem#2 are plotted at two different time 
(a) t = 2; (b) t = 4. 


27 





























































(a) 



(b) 

Figure 9: Magnified solution for the exponential wave at t = 4. (a) Extreme; (b) Bottom. 
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(b) 

Figure 10: Magnified solution for the square waves at t = 4. (a) Extreme; (b) Bottom. 
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(a) 



(b) 

Figure 11: The predieted results for linear adveetion problem#2 are plotted at t = 4. 
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(a) 



(b) 

Figure 12: Comparison of the results using the present hybrid strategy and the present hybrid 
strategy (a) rc = 0.1; (b) = 10.0. 


31 












































(a) (b) 



(c) (d) 


Figure 13: (a) The predicted results for linear advection problem#3 are plotted at t = 20; 
(b)Magnified solution between —1.4 < jc < —0.4; (c) Magnified solution between —0.7 < .r < 
0.1; (d)Magnified solution between 0 < jc < 0.7. 
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(b) 

Figure 14: The predieted results for non-linear adveetion problem are plotted at two different 
time, (a) t = 0.3; (b) t = 0.35. 
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(C) (d) 

Figure 15: The predieted results for the Shu-Osher problem are plotted at four different time, 
(a) t = 0.45; (b) t = 0.9; (e) t = 1.35; (d) t = 1.8. 
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Figure 16: Comparison of the results by ODRPCCD5 and WEN05 sehemes for the vortex flow 
problem eomputed in 100 x 100 grids, (a) t = 1.5; (b) t = 2.5; (e) t = 4; (d) t = 5. 
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Figure 17: Comparison of the results by ODRPCCD5 and WEN05 sehemes at t = 2.5. 
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Figure 18: The predicted results for the Zalesak’s problem, (a) 100 x 100 grids; (b) 200 x 200 
grids. 
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